Our goal will be to produce an equivalence table relating year 2010 Bay Area Census Tracts to MTC’s Transportation Analysis Zones (TAZ).
Mapview is only required if you want to render maps.
dplyr, readr, and sf are required by the analysis but sf could probably be removed if necessary.
library(sf)
library(dplyr)
library(readr)
library(mapview)
Please see the data processing doc (data_prep.Rmd) for background on data sources and calculated variables.
#setwd("~/Box/")
#setwd("~/Box/DataViz Projects/")
setwd("~/Box/DataViz Projects/Data Analysis and Visualization/census_examples/taz_tract/data")
blocks <- st_read("blocks.gpkg")
tracts10 <- st_read("tracts.gpkg")
tracts00 <- st_read("tracts00.gpkg")
taz_df <- st_read("taz1454.gpkg")
intersection_df <- st_read("intersection_df.gpkg")
intersection_df$probably_a_sliver <- as.logical(intersection_df$probably_a_sliver)
intersection_df$definitely_a_sliver <- as.logical(intersection_df$definitely_a_sliver)
Tract_zone_2000 <- read_csv("Tract_zone_2000.csv")
Tract_zone_2000 <- rename(Tract_zone_2000, tract = Tract)
year_2000_intersection_df <- st_read("intersection_df_2000.gpkg")
From review of the year 2000 lookup table in other scripts, we’ve defined the goals of the table:
For any given Tract be able to look up either: (a) the TAZ that fully and completely circumscribes it, (b) that mostly circumscribes it, or finally, (c) the set of TAZ’s that are within it.
We have a rough idea of how many intersections we should expect to fall in each category by looking at the year 2000 CSV.
tract_zone_2000_dense <- reshape2::melt(Tract_zone_2000, id.vars = "tract", na.rm=TRUE)
number_of_equivalencies <- group_by(tract_zone_2000_dense, variable) %>%
summarize(count = n())
total <- sum(number_of_equivalencies$count)
number_of_equivalencies$percent <- round(number_of_equivalencies$count/total, digits=3)
print(number_of_equivalencies)
## # A tibble: 5 x 3
## variable count percent
## <fct> <int> <dbl>
## 1 rtaz1 1405 0.966
## 2 rtaz2 30 0.0210
## 3 rtaz3 11 0.00800
## 4 rtaz4 5 0.00300
## 5 rtaz5 3 0.00200
rm(number_of_equivalencies)
We start by build the equivalence for new tracts only, assuming the 2000 table is correct.
‘rtaz’ is a column header from the year 2000 lookup table. If its not on the tract, then the tract hasn’t already been mapped to a TAZ.
new_tract_ids <- tracts10[is.na(tracts10$rtaz1),]$tract
So we have 387 new tracts in 2010. By numbers this might seem surprising. There are 1.4k tracts in 2000 and 1.5k in 2010. Definitely not close to 400 new tracts by number. Lets have a look at the tracts that are new for 2010, and those from 2000 that aren’t in 2010 to understand what changed.
tracts00_ids <- tracts00[!tracts00$tract %in% tracts10$tract,]$tract
tracts10_ids <- tracts10[!tracts10$tract %in% tracts00$tract,]$tract
mapview(tracts10[tracts10$tract %in% tracts10_ids,], col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(tracts00[tracts00$tract %in% tracts00_ids,], col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light'))
Looking at this its clear that some tracts were split, and others may have had other changes. A detailed list of changes made is available here:
https://www.census.gov/geo/reference/bndrychange/changenotedisplay.php
For our purposes, we’ll just update the 387 tracts with new IDs, making a note to drop the tracts no longer in 2010 by ID when we construct a new table.
tracts_to_drop <- tracts00_ids
So first we will drop the “probably a sliver” spatial relationships. This includes “definitely a sliver” relationships. These are not relationships that fall in any of the above categories.
intersection_df_new_tracts <- intersection_df[!intersection_df$probably_a_sliver==TRUE,]
intersection_df_new_tracts <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% new_tract_ids,]
dim(intersection_df_new_tracts)
## [1] 509 14
That leaves us with about 500 relationships for for 300 tracts.
In practice, there is at least one error in the 2000 lookup table that we will want to update.
The second taz (rtaz) for tract 500300 is incorrect it should. be mapped to TAZ 539 but was mapped to TAZ 529. More details at the bottom here: https://bayareametro.github.io/Data-And-Visualization-Projects/census_examples/taz_tract/reverse_engineer_tract_zone_2000_method.html
We’ll fix that here since we will use the 2000 data as a basis for 2010.
print(Tract_zone_2000[Tract_zone_2000$tract=="500300",'rtaz2'])
## # A tibble: 1 x 1
## rtaz2
## <int>
## 1 529
Tract_zone_2000[Tract_zone_2000$tract=="500300",'rtaz2']=539
print(Tract_zone_2000[Tract_zone_2000$tract=="500300",'rtaz2'])
## # A tibble: 1 x 1
## rtaz2
## <dbl>
## 1 539.
Lets first identify all tracts in sets (a) and (b)
These tracts should have an area of intersection with a TAZ which is at least as large (within error) as the TAZ.
We say “within error” because the topology errors here make full identity unlikely.
Ideally, the Tract’s area of intersection with the TAZ would be identical to the TAZ’s area (or greater). But because of topological errors, this rarely occurs. So we have to decide what cutoff we want to say is “within error.”
Below we set it to conditions where the Tract’s area of intersection with the TAZ is at least 70% of the area of the TAZ.
We can compare this to the ratio of single-taz intersections (types (a) and (b)) in the year 2000 data.
intersection_df_new_tracts$intersection_area_over_tract_area <- as.numeric(intersection_df_new_tracts$intersection_area)/as.numeric(intersection_df_new_tracts$tract_area)
summary(intersection_df_new_tracts$intersection_area_over_tract_area)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0001366 0.3553154 0.9610558 0.7221450 0.9862458 1.0000000
print(table(intersection_df_new_tracts$intersection_area_over_tract_area>.70))
##
## FALSE TRUE
## 146 363
If we base our decision on what on 2000 data then we should expect about 96% of the Tracts to map to only 1 TAZ, which again would be the relationship conditions defined as (a) and (b).
At 70% of intersection we get 363/387 or 93% of the intersections in the new Tracts being single-taz intersections.
This seems about right, and a bit conservative.
Interestingly, setting the intersection_area_over_tract_area ratio to 50% gets us to almost exactly 96%, which is what the year 2000 data had. But lets be conservative and go with 70%
Lets call those intersections acceptable, and put them in their own data frame.
intersection_type_a_b <- intersection_df_new_tracts[intersection_df_new_tracts$intersection_area_over_tract_area>.70,]
tracts_left_id <- intersection_df_new_tracts[!intersection_df_new_tracts$tract %in% intersection_type_a_b$tract,]$tract
length(tracts_left_id)
## [1] 35
We’ll leave that data frame (intersection_type_a_b) aside and add them to the main table later with the type (c) intersections.
That leaves us just a handful of tracts left. We’ll put the spatial intersections for these types (c) into their own data frame.
intersection_tracts_left <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tracts_left_id,]
dim(intersection_tracts_left)
## [1] 35 15
rm(intersection_tracts_left)
Lets have a look at them:
tract_ids <- tracts_left_id
taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz
m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]
mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))
Below we’ll review a few of these intersections in more detail in order to determine which we should keep.
Others are examples of where the TAZ is larger than the Tract. For example:
cross_tract <- "505009"
tract_ids <- cross_tract
taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz
m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]
mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))
This is a type (b) and type (c) coverage. Some of the remaining tracts are like this. From here on out we’ll just call the intersections valid or not.
In some cases the Tract seems to be split by a TAZ.
cross_tract <- "403502"
tract_ids <- cross_tract
#m_blocks <- blocks[blocks$TRACTCE10 %in% tract_ids,]
taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz
m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]
mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))
We might want to develop rules for this kind of TAZ/Tract relationship in the future. For now, we will just say its equivalent to the 1 TAZ that majority covers it (972).
intersection_type_a_b <- rbind(intersection_type_a_b,
intersection_df_new_tracts[intersection_df_new_tracts$tract=="403502" &
intersection_df_new_tracts$taz==972,])
Valid intersections include those in tracts: 511705, 614000, 061500, 061100, 131100, 410500, 313206, 061200,601902, 601901
990100 is not valid, because it is a park/open space preserve. 980401 is not valid because its the Farallon islands.
Lets add the valid intersections to the table.
valid_intersections <- c("505009","613900","061500","511705", "614000",
"061500", "061100", "131100", "410500",
"313206","061200","601902", "601901")
intersection_type_c <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% valid_intersections,]
intersection_type_c <- intersection_type_c
st_geometry(intersection_type_c) <- NULL
intersection_type_c <- dplyr::select(intersection_type_c,taz,tract)
intersection_type_c <- intersection_type_c[,c('tract','taz')]
intersection_type_c$header_string <- 'rtaz'
intersection_type_c$num <- ave(intersection_type_c[['taz']],
intersection_type_c[['tract']],
FUN = seq_along)
intersection_type_c_sparse <- intersection_type_c %>%
tidyr::unite("header_string",
header_string,
num) %>%
tidyr::spread(header_string, taz)
#rm(intersection_type_c)
When we build the sparse table for multiple intersections, it becomes apparent that we will have to add a column to the lookup table. In 2000, there were 5 columns, but here we’ve found that (for 2010) there are 6 valid intersections for 1 tract. That tract is in downtown SF, and a map of it can be found in the Appendix.
Lets take a look at what we’re adding to the updated 2000 table before we do it. We may need to reshape the columns.
st_geometry(intersection_type_a_b) <- NULL
intersection_type_a_b <- intersection_type_a_b[,c('tract','taz')]
knitr::kable(head(intersection_type_c_sparse))
| tract | rtaz_1 | rtaz_2 | rtaz_3 | rtaz_4 | rtaz_5 | rtaz_6 |
|---|---|---|---|---|---|---|
| 061100 | 25 | 24 | NA | NA | NA | NA |
| 061200 | 143 | 141 | NA | NA | NA | NA |
| 061500 | 12 | 13 | 15 | 14 | 16 | 17 |
| 131100 | 1417 | 1453 | 1454 | NA | NA | NA |
| 313206 | 1202 | 1201 | NA | NA | NA | NA |
| 410500 | 985 | 984 | NA | NA | NA | NA |
knitr::kable(head(intersection_type_a_b))
| tract | taz | |
|---|---|---|
| 77 | 503337 | 651 |
| 80 | 503336 | 651 |
| 84 | 503330 | 653 |
| 85 | 503329 | 653 |
| 338 | 508003 | 479 |
| 339 | 508004 | 479 |
knitr::kable(head(Tract_zone_2000))
| tract | rtaz1 | rtaz2 | rtaz3 | rtaz4 | rtaz5 |
|---|---|---|---|---|---|
| 010100 | 41 | NA | NA | NA | NA |
| 010200 | 40 | NA | NA | NA | NA |
| 010300 | 39 | NA | NA | NA | NA |
| 010400 | 38 | NA | NA | NA | NA |
| 010500 | 22 | 23 | NA | NA | NA |
| 010600 | 37 | NA | NA | NA | NA |
So we need to reshape the type_a_b table and add the new column to the 2000 table.
And rename columns to match (underscores, etc)
Tract_zone_2010 <- Tract_zone_2000
intersection_type_a_b <- dplyr::rename(intersection_type_a_b, rtaz1 = taz)
intersection_type_a_b$rtaz2 <- as.integer(NA)
intersection_type_a_b$rtaz3 <- as.integer(NA)
intersection_type_a_b$rtaz4 <- as.integer(NA)
intersection_type_a_b$rtaz5 <- as.integer(NA)
intersection_type_a_b$rtaz6 <- as.integer(NA)
Tract_zone_2010$rtaz6 <- as.integer(NA)
names(intersection_type_a_b) <- names(Tract_zone_2010)
names(intersection_type_c_sparse) <- names(Tract_zone_2010)
Now we can bind them all together. Lets have a look at them and then do so.
knitr::kable(head(Tract_zone_2010))
| tract | rtaz1 | rtaz2 | rtaz3 | rtaz4 | rtaz5 | rtaz6 |
|---|---|---|---|---|---|---|
| 010100 | 41 | NA | NA | NA | NA | NA |
| 010200 | 40 | NA | NA | NA | NA | NA |
| 010300 | 39 | NA | NA | NA | NA | NA |
| 010400 | 38 | NA | NA | NA | NA | NA |
| 010500 | 22 | 23 | NA | NA | NA | NA |
| 010600 | 37 | NA | NA | NA | NA | NA |
knitr::kable(head(intersection_type_a_b))
| tract | rtaz1 | rtaz2 | rtaz3 | rtaz4 | rtaz5 | rtaz6 | |
|---|---|---|---|---|---|---|---|
| 77 | 503337 | 651 | NA | NA | NA | NA | NA |
| 80 | 503336 | 651 | NA | NA | NA | NA | NA |
| 84 | 503330 | 653 | NA | NA | NA | NA | NA |
| 85 | 503329 | 653 | NA | NA | NA | NA | NA |
| 338 | 508003 | 479 | NA | NA | NA | NA | NA |
| 339 | 508004 | 479 | NA | NA | NA | NA | NA |
knitr::kable(head(intersection_type_c_sparse))
| tract | rtaz1 | rtaz2 | rtaz3 | rtaz4 | rtaz5 | rtaz6 |
|---|---|---|---|---|---|---|
| 061100 | 25 | 24 | NA | NA | NA | NA |
| 061200 | 143 | 141 | NA | NA | NA | NA |
| 061500 | 12 | 13 | 15 | 14 | 16 | 17 |
| 131100 | 1417 | 1453 | 1454 | NA | NA | NA |
| 313206 | 1202 | 1201 | NA | NA | NA | NA |
| 410500 | 985 | 984 | NA | NA | NA | NA |
Tract_zone_2010 <- rbind(Tract_zone_2010,intersection_type_a_b,intersection_type_c_sparse)
dim(Tract_zone_2010)
## [1] 1781 7
#rm(intersection_type_a_b)
#rm(intersection_type_c_sparse)
We have too many tracts now. Thats because we still have tracts that were changed from 2000 to 2010 in the table that we have not dropped.
Drop the tracts updated for 2010 (by id)
Tract_zone_2010 <- Tract_zone_2010[!Tract_zone_2010$tract %in% tracts_to_drop,]
dim(Tract_zone_2010)
## [1] 1581 7
dim(tracts10)
## [1] 1592 10
We still have a few missing.
Map those tracts that are missing from the lookup.
missing_tracts10 <- tracts10[!tracts10$tract %in% Tract_zone_2010$tract,]
These should just be water tracts. We will load tracts not cut by water and water here to double check that.
options(tigris_use_cache=TRUE)
library(tigris)
counties1<-c("01","13","41","55","75","81","85","95","97")
tracts_not_cut <- tigris::tracts("CA", counties1, class="sf", year=2010)
detach("package:tigris", unload=TRUE)
tracts_not_cut <- dplyr::select(tracts_not_cut,TRACTCE10)
tracts_not_cut <- dplyr::rename(tracts_not_cut,tract = TRACTCE10)
tracts_not_cut <- st_transform(tracts_not_cut, crs=26910)
tracts_not_cut$tract_total_area <- st_area(tracts_not_cut)
bay_water <- st_read("https://geo.nyu.edu/download/file/stanford-mb777jk0330-geojson.json")
## Reading layer `stanford-mb777jk0330-geojson' from data source `https://geo.nyu.edu/download/file/stanford-mb777jk0330-geojson.json' using driver `GeoJSON'
## Simple feature collection with 367 features and 6 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -123.6336 ymin: 36.92222 xmax: -121.2531 ymax: 38.85059
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
bay_water <- bay_water[st_is_valid(bay_water),]
bay_water <- st_transform(bay_water, crs=26910)
df1 <- st_intersection(tracts_not_cut, bay_water)
df1$intersection_area <- st_area(df1)
mapview(df1[df1$tract %in% missing_tracts10$tract,])
rm(bay_water)
rm(tracts_not_cut)
rm(df1)
These tracts are all over water.
We can put them into the table for clarity
water_tracts10 <- data.frame(tract=c(unique(as.character(missing_tracts10$tract))))
water_tracts10$rtaz1 <- as.integer(NA)
water_tracts10$rtaz2 <- as.integer(NA)
water_tracts10$rtaz3 <- as.integer(NA)
water_tracts10$rtaz4 <- as.integer(NA)
water_tracts10$rtaz5 <- as.integer(NA)
water_tracts10$rtaz6 <- as.integer(NA)
knitr::kable(head(water_tracts10))
| tract | rtaz1 | rtaz2 | rtaz3 | rtaz4 | rtaz5 | rtaz6 |
|---|---|---|---|---|---|---|
| 990100 | NA | NA | NA | NA | NA | NA |
| 990000 | NA | NA | NA | NA | NA | NA |
| 980401 | NA | NA | NA | NA | NA | NA |
Tract_zone_2010 <- rbind(Tract_zone_2010,water_tracts10)
dim(Tract_zone_2010)[1]==length(unique(tracts10$tract))
## [1] TRUE
rm(water_tracts10)
It has the right number of rows
tract_zone_2010_dense <- reshape2::melt(Tract_zone_2010, id.vars = "tract", na.rm=TRUE)
tazs <- unique(tract_zone_2010_dense$value)
table(taz_df$taz %in% tazs)
##
## FALSE TRUE
## 8 1446
There are 8 TAZ’s not in the table. Were they in the intersection types a, b or c we identified?
missing_tazs <- taz_df[!taz_df$taz %in% tazs,]$taz
table(missing_tazs %in% intersection_type_a_b$rtaz1)
##
## FALSE
## 8
table(missing_tazs %in% intersection_type_c$taz)
##
## FALSE
## 8
They were not.
Are they not new tracts?
table(intersection_df[intersection_df$taz %in% missing_tazs,]$tract %in% tracts00$tract)
##
## FALSE TRUE
## 11 13
They were not because they were in the year 2000 tract table.
tract_zone_2000_dense[tract_zone_2000_dense$taz %in% missing_tazs,]
## [1] tract variable value
## <0 rows> (or 0-length row.names)
They were not in the 2000 equivalency table.
Where are they?
tract_overlap_ids <- intersection_df[intersection_df$taz %in% missing_tazs,]$tract
m_tracts10 <- tracts10[tracts10$tract %in% tract_overlap_ids,]
m_taz <- taz_df[taz_df$taz %in% missing_tazs,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_overlap_ids,]
mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))
These seem like intersections that should be included based on our definitions. All of these would fall into both types (b) and (c) intersections.
So we’ll add them to the table
additions_to_2000_lookup <- m_intersection[,c('tract','taz')]
st_geometry(additions_to_2000_lookup) <- NULL
Tract_zone_2010_temp <- Tract_zone_2010[!Tract_zone_2010$tract %in% additions_to_2000_lookup$tract,]
additions_to_2000_lookup$header_string <- 'rtaz'
additions_to_2000_lookup$num <- ave(additions_to_2000_lookup[['taz']],
additions_to_2000_lookup[['tract']],
FUN = seq_along)
additions_to_2000_lookup_sparse <- additions_to_2000_lookup %>%
tidyr::unite("header_string",
header_string,
sep="",
num) %>%
tidyr::spread(header_string, taz)
additions_to_2000_lookup_sparse$rtaz6 <- as.integer(NA)
Tract_zone_2010_temp <- rbind(Tract_zone_2010_temp,additions_to_2000_lookup_sparse)
dim(Tract_zone_2010_temp)
## [1] 1584 7
Check again for missing TAZ’s:
tract_zone_2010_dense <- reshape2::melt(Tract_zone_2010_temp, id.vars = "tract", na.rm=TRUE)
tazs <- unique(tract_zone_2010_dense$value)
table(taz_df$taz %in% tazs)
##
## TRUE
## 1454
No missing TAZ’s.
missing_tracts <- tracts10[!tracts10$tract %in% Tract_zone_2010$tract,]
dim(missing_tracts)
## [1] 0 10
No missing tracts.
In sum, we updated 1 typo in the 2000 table, 10-20 intersections that should have been included, and updated 387 new/split tracts with their meaningful spatial intersections as defined at the outset.
setwd("~/Documents/Projects/mtc/Data-And-Visualization-Projects/census_examples/taz_tract")
write_csv(Tract_zone_2010,'Tract_zone_2010.csv')
The relationships, in terms of areas, are rarely “equivalent”. Furthermore, the kind of equivalence in area categeorically varies across rows, so the way that one tract is “equivalent” to a TAZ is not the same as the way another tract is “equivalent” to a TAZ.
For example, tracts can be smaller than a TAZ and TAZ’s can be smaller than a tract. In some cases, a tract is called “equivalent” to a tract when it is 1/20th the size of that TAZ. In other cases, a TAZ is called “equivalent” to a Tract when it is 1/20th the size of that Tract.
It seems that for modeling, there is an understanding of “equivalence” which may be worth describing in another document.
All TAZ’s can be made spatially equivalent to all Blocks, but not all TAZ’s can be made spatially equivalent to Tracts
It may be useful to think about these equivalencies in terms of smaller spatial units. When we compare a few of the TAZ’s to blocks it is apparent that they were built on blocks, not tracts.
Therefore, the more accurate “equivalence” for all TAZ’s is described at the block (or perhaps block group) level, not at the tract.
While it is true that the majority of TAZ’s can be built from Tracts, many cannot.
For this small set, there is a loss of accuracy in aggregating them up a geographic level.
On the other hand, there would be no loss in accuracy in disaggregating the measures at the Tracts level down to blocks, as long as the error and disaggregation is represented properly.
Two example tracts follow. There are other instances in which a TAZ to Block lookup would make more sense, especially when tracts are redrawn.
The first example is also the Tract that required adding a new column to the table.
block_tract1 <- "061500"
tract_ids <- block_tract1
m_blocks <- blocks[blocks$TRACTCE10 %in% tract_ids,]
taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz
m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]
mapview(m_tracts10, color="red", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))
mapview(m_taz, color="blue", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))
mapview(m_intersection,color="green", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))
mapview(m_blocks, color="magenta", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))
mapview(m_tracts10, color="red", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light')) +
mapview(m_blocks, color="magenta", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light')) +
mapview(m_taz, color="blue", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light'))
If TAZ boundaries can be built from Census Blocks we might need not have to draw geometries at all.
Instead, Travel Model geographies could reference census defined areas at the lowest level of aggregation.
Then we just rely on Census to maintain the geometries with no loss of information/accuracy.